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ABSTRACT 

Two new equivalent linearization implementations for geometrieally nonlinear random 
vibrations are presented. Both implementations are based upon a novel approaeh for evaluating 
the nonlinear stiffness within commercial finite element codes and are suitable for use with any 
finite element code having geometrically nonlinear static analysis capabilities. The formulation 
includes a traditional force-error minimization approach and a relatively new version of a 
potential energy-error minimization approach, which has been generalized for multiple degree- 
of-freedom systems. Results for a simply supported plate under random acoustic excitation are 
presented and comparisons of the displacement root-mean-square values and power spectral 
densities are made with results from a nonlinear time domain numerical simulation. 

1. INTRODUCTION 

Current efforts to extend the performance and flight envelope of high-speed aerospace vehicles 
have resulted in structures which may respond to the imposed dynamic loads in a geometrically 
nonlinear (large deflection) random fashion. What differentiates the geometrically nonlinear 
random response considered in this paper from a linear response is the presence of bending- 
membrane coupling, which gives rise to membrane stretching (in-plane stresses) in the former. 
This coupling has the effect of stiffening the structure and reducing the dynamic response 
relative to that of the linear system. Linear analyses do not account for this effect and 
consequently may significantly over-predict the response, leading to grossly conservative 
designs. Without practical design tools capable of capturing the nonlinear dynamics, further 
improvements in vehicle performance and system design will be hampered. 

Methods currently used to predict geometrically nonlinear random response include perturbation, 
Fokker-Plank-Kolmogorov (F-P-K), numerical simulation and stochastic linearization 
teehniques. All have various limitations. Perturbation techniques are limited to weak geometrie 
nonlinearities. The F-P-K approach [1, 2] yields exact solutions, but can only be applied to 
simple meehanieal systems. Numerical simulation techniques using numerieal integration 
provide time histories of the response from which statistics of the random response may be 
ealeulated. This, however, comes at a high computational expense due to the long time reeords 
or high number of ensemble averages required to get high quality random response statisties. 
Statistieal linearization methods, for example equivalent linearization (EL) [[2-7]), have seen the 
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broadest application because of their ability to accurately capture the response statistics over a 
wide range of response levels while maintaining relatively light computational burden. In EL 
methods, an equivalent linear system is sought such that the difference between the equivalent 
linear and nonlinear systems is minimized. Traditionally, the difference between the nonlinear 
force and the product of the equivalent linear stiffness and displacement response vector is 
minimized. 

Since the analysis of complicated structural geometries is required, a finite element-based EL 
method is considered most appropriate. In the past, implementations of EL using finite element 
analysis have been limited to special purpose codes. This is largely due to the inaccessibility of 
the nonlinear element formulation in commercial finite element applications. For example, the 
nonlinear stiffness of a rectangular plate finite element was found in [8] after a lengthy analytical 
derivation, which utilized information about the element shape functions. 

The first known EL implementation in a general-purpose finite element code [9] was developed 
for use in MSC.NASTRAN [10] version 67. In that implementation, called “Super Element 
Modal Equivalent Linear Random Response” or SEMELRR, the equivalent linear stiffness was 
obtained as the sum of the linear stiffness and three times the differential stiffness. While this 
form was convenient to implement, it was found to over-predict the degree of nonlinearity and 
produce non-conservative results. Over-prediction of nonlinearity can produce the undesirable 
result of structural designs incapable of withstanding the applied loads in an acceptable fashion. 

This paper describes a novel approach to accurately determine the nonlinear stiffness without 
using a description of the element shape functions [11, 12]. The approach solves a series of 
inverse linear and nonlinear static problems to evaluate the nonlinear force, from whieh the 
nonlinear stiffness can be determined. While it is applicable to any commercial finite element 
program having a nonlinear static analysis capability, MSC.NASTRAN was selected due to its 
widespread use in the aerospace industry. The stiffness evaluation technique was first validated 
for a beam strueture under a special loading condition [12]. More recently, it was validated for a 
elamped-elamped beam under random inertial loading through comparison with a finite element- 
based numerieal simulation analysis in physical coordinates [13]. 
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In addition, this paper describes the generalization of a relatively new potential energy-error 
minimization version of the EL method [6, 7] to a multiple degree-of-freedom case. This new 
minimization approach was previously validated using the F-P-K method for a few special cases 
including a Duffing oscillator and a two degree-of-freedom system [11] and for a beam structure 
[ 12 ]. 

To further validate the nonlinear stiffness evaluation approach and its use in an EL analysis using 
force-error and potential energy-error minimization, a simply supported plate subjected to 
random acoustic loading is considered. This case is also closely related to the intended 
application. Results are compared with a finite element-based numerical simulation analysis in 
modal coordinates, as described in [14]. The numerical simulation analysis in modal coordinates 
was selected for this purpose because the solution in physical coordinates was intractable due to 
the large system size. The use of a modal approach is considered acceptable for the class of 
problems exhibiting bending-membrane coupling, as in the simply supported plate. 

2. NONLINEAR STIFFNESS EVALUATION 

The equations of motion of a multiple degree-of-lfeedom, viscously damped geometrically 
nonlinear system can be written in the form: 

MX{t) + cx{t) + KX{t) + r(Ar(0) = F{t) (i) 

where M, C, K are the mass, damping, and stiffness matrices, X is the displacement response 
vector and F is the force excitation vector, respectively. The nonlinear force term r(A) is a 
vector function, which generally includes second and third order terms in A. 

Solution to equation (1) via any method requires knowledge of the system matrices. In the 
eontext of a eommercial finite element program, M, C, and K are generally available. In the EL 
analyses to follow, the nonlinear stiffness is required. The nonlinear stiffness is related toE , 
whieh is typieally not available within a commercial finite element program. Therefore, a means 
of numerieally evaluating E was developed, as is next described, for the determination of the 
nonlinear stiffness. 

A set of eoupled modal equations with reduced degrees-of-lfeedom is first obtained by applying 
the modal eoordinate transformation 

X = ^q ( 2 ) 
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to equation (1), where D is generally a subset {L < N) of the linear eigenveetors , q is the veetor 
of modal eoordinates, and the time dependence is implied. This coupled set is expressed as 

MqWCqWkqW':i{q^,q^,...,qi}W F (3) 

where 

MDD^MD DT/J 

CDD^CD W\lC,co^\ 

D W\col\ 

Y D D 
FWW^F 

the components of q, and O)^ are the undamped natural frequencies. The 
components of the nonlinear force vector may be written in the form 

!. L L L L 

Yj Y ^jk^j^k Q X S Z ^jki<ij^k<li r D 1, 2, . . . , L . (4) 

,/Dl kUj /Dl kOj /Di 

where a'j^. and are nonlinear stiffness coefficients with j = \,1, L, k =j,j+\, ...,L and / 
= k, k+\, L. This particular form of y facilitates the subsequent solution of the equivalent 
linear system. Its evaluation entails solving for the coefficients aj,^ and bj^, using a new 

procedure developed for use with finite element programs having a nonlinear static solution 
capability. 

The procedure is based on the restoration of nodal applied forces by prescribing nodal 
displacements to both linear and nonlinear static solutions. The total nodal force F^ may be 
written in physical coordinates as 

f,Df,Df„D«a:,DD(a:j (5) 

where is a prescribed physical nodal displacement vector, and F^^ and the linear and 

nonlinear contributions to the total nodal force. is first obtained by prescribing in the 
linear static solution. F^ is then obtained by prescribing X^ in the nonlinear static solution 
which includes both linear and nonlinear contributions. Finally, the nonlinear contribution F^^^ 
is obtained by subtracting F^ from F^, , or 
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( 6 ) 


To illustrate the teehnique, one can begin by prescribing the displacement fields 

X, D 

D 

The nonlinear nodal force contributions Fj^j^ are determined using ( 6 ) after solving the linear and 
nonlinear static solutions. These may be written in modal coordinates as 

Fnl, D D ^Fnl, D D D ^[n^<h<h<h 

Fnl, D D ^Fnl, D D D 

where the sought stiffness coefficients {a[^\ and [ 6 ^,] are column vectors of length L. Note that 
the other nonlinear terms do not appear in (7) since qj DO for j Since is a known scalar, 
the coefficients [a[J and for r D 1,2, ,L can be determined from the resulting system (7) 
of 2 x 1 , linear equations. The remaining coefficients and (7 D 2,3 ,L) can be 
determined in an analogous manner. 

A similar technique can be employed to determine stiffness coefficients with two unequal lower 
indices, e.g., and [ 6 , 22 ]- Coefficients of this type appear only if the number of 

retained eigenvectors is greater than or equal to two (L > 2). Prescribing the displacement fields 

D D 4 >,g, D %^2 
D D 4 »,g, - 4^^2 

results in the following equations 


5 



Fnl, D D D 4^2) D §«n§9l^l D §^ni§^l^l^l D §«22§^2^2 D §422§^2^2^2 

D §«i' 2§M2 D §^n2§^1^1^2 D §422§^1^2^2 

Fnl, D D ^D(D4^1 D 4^2) D §4i§?i9i D I^Uli^l^l^l D §«22§^2^2 ^ §44 | ^2 ^2 ^2 

D §«r2§M2 D §^n2§^1^1^2 D §422§^1^2^2 
FnL, D D ^0(0 4^1 D 4^2) D §<1§?1^1 D §^ni§^l^l^l D §«22§^2^2 D §422§^2^2^2 

D §«i'2§^i 92 D §^n2§^1^1^2 D §422§^1^2^2 

Summing the first two of equations (8) results in 

Fnu D 4vi, D 2 §«;, D 2 ^q^q^ D 2 ^q,q^ 


from whieh the eoeffieients [ajj] may be determined, since [a[]] and [ 022 ] were previously 
determined. Then, from the first and third of equations (8), the coefficients [h,''j 2 ] and [^[22 ] may 
be determined from the 2 D Z, system of equations. In this manner, all coefficients of the type 
Wjj^'\ and for j,k D 1,2, ,L may be determined. 

For cases when the number of retained eigenvectors is greater than or equal to three {L D 3), 
coefficients with three unequal lower indices, e.g., [44 ]’ t>e determined by prescribing the 
displacement field 

□ □4^1 D 4^2 D 4% • 


The resulting equation 

Fnl D D (04^1 D 4^2 D 4 ^ 3 ) 

D §«n§^l^l D §«22§^2^2 D §^33§^393 ^ §«l''2 §9l^2 ^ ^ 

(9) 

D §4li§^l^l^l D §^24§^2^292 D §433 §93^3^3 ^ §^H2 §^l9l^2 ^ §44 §92?29l 
^ §4 i3§?1^1^3 ^ §431§?3?3^1 ^ §^223 § 92^2^3 ^ §^332 §93^392 ^ §423§^1^2^3 

contains one column of unknown coefficients [441 • coefficients of type [hyt/Ki ^ kW 1) 
can be found in this manner. 


Having the modal equations of motion (3) formulated, their solution can be undertaken through a 
variety of techniques. For the reasons previously discussed, application of EL is considered. 
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3. EQUIVALENT LINEARIZATION APPROACH 

An approximate solution to (1) can be achieved by formation of an equivalent linear system: 

MX{t) D CX{t) D (A D AJA(0 D F{t) (10) 

where is the equivalent linear stiffness matrix. While it is possible to perform an EL analysis 

in the physical degrees-of-freedom, it is desirable to recast the problem in modal coordinates to 
simplify the problem. The equivalent linear analog of equation (3) may be found by applying the 
modal transformation (2) to (10): 

MqWCqW^WK^^qW F (11) 

where the fully populated modal equivalent stiffness matrix is given by 

A.D DXa • 

Two EL approaches are considered. One is based on minimization of the error in the nonlinear 
force vector and the other minimizes the error in potential (strain) energy. 

3.1. FORCE-ERROR MINIMIZATION APPROACH 

The traditional (force-error minimization) method of EL seeks to minimize the difference 
between the nonlinear force and the product of the equivalent linear stiffness and displacement 
response vector. Since the error is a random function of time, the required condition is that the 
expectation of the mean square error be a minimum [3], i.e. 

error D E p (X) D K^XfiQ (X) D A,A)§ ^ min (12) 

where A[...] represents the expectation operator. Equation (12) will be satisfied if 


d(error) 


D 0 


hyDU, ,N 


where are the elements of matrix , and N is the number of physical degrees of freedom. 


In this study, consideration is limited to the case of Gaussian, zero-mean excitation and response 
to simplily the solution. With these assumptions and omitting intermediate derivations, the final 
form for the equivalent linear stiffness matrix in physical coordinates becomes (see for example 
[4], [5]): 
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and in modal coordinates becomes 


K^uh 

D 


ODD DP 


with D as previously defined in (4). 


(14) 


3.2. POTENTIAL ENERGY-ERROR MINIMIZATION APPROACH 

An alternative EL approach based on potential (strain) energy-error minimization was proposed 
in [6, 7]. Analysis in these works was limited to single degree-of-lfeedom systems and 
simplified multiple degree-of-freedom systems. In this study, that approach is rigorously 
generalized for the case of coupled multiple degree-of-freedom systems. One can begin with an 
expression for the error in potential energy 

error D £g(£/(;v)0+x'-jir,;v(g (is) 


where U(X) is the potential energy of the original (nonlinear) system. A condition of 
minimized error requires that 


D 




-g(f/(A)DiA^IST,A)' 


D 0 


h7Ql,2, ,N. 


(16) 


Omitting intermediate derivations, one obtains the following system of linear equations with 
respect to unknown elements of matrix : 


D D K^..E^x.x.x^x, 


;Dl ./Dl 


D 2E[x^XiU{X)] 


k,/Dl,2, ,N. 


(17) 


For example, equation (17) would have the following form for a two degree-of-freedom system 



D fill;! 

E^lxl'^ Ey,xl 

file, 41 e^A 




■1 

|£|»rfC/(X)|[ 

0' 

Ki2[ 

]D2[ 

]£Dt,x,y(Ar)[I 

Kn[ 

‘f[ic,xjiJ(X)[- 


i 

]£|x;£/(X)|[ 


( 18 ) 


This system of equations is under-determined (the second and third rows are identieal), so an 
additional equation is required to solve this system. The additional equation(s) can be provided 
by imposing symmetry of the matrix K^, i.e., K^y D K^y. 


Note that the formulation for the multiple degree-of-freedom system in reference [7] differs from 
the above formulation. The multiple degree-of-freedom case in [7] is treated by introduction of 
only diagonal terms of the matrix K^, whereas in this study the complete matrix is 
formulated and is thereby advantageous in cases which exhibit modal coupling. 


In modal coordinates, the equivalent linear stiffness matrix is related to nonlinear stiffness terms 
in a more complicated manner than in the force-error minimization approach. It is known that 
the nonlinear elastic force terms satisfy the following 






wu 


Wq, 


rDl,2, ,L 


(19) 


where U is the potential energy generated by nonlinear terms only. The potential energy of the 
system may be written in the following form 


[/ □ D D D D dsjkiqsqjqkqi (20) 

.?Dl /D.V kU J lUk 

whieh upon substitution into (19) yields 


Uq^ 


L L L L 


D D D D dy^q^ .q^q, 

.yDl /D.y kW j iWk 


( 21 ) 


The eoeffieients dy^., are related to the nonlinear stiffness coefficients through equation (4) as: 


D D D D D D D 

,/Dl kWj ./Dl kWj ;Di 




L L L L A 

D D D D dsjkiqxqjqkqi 

.yDl jUs kU J lUk J 


( 22 ) 
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In the following, a zero-mean response is assumed, i.e., E{q\ DO. As will be subsequently 
shown, this has the effeet of dropping out the quadratic terms . From (22), the coefficients 
can then be written as: 


^/2 

sW j 


5 D j W k 


5 D 7 D ^ D / 


otherwise. 


(23) 


Note that if the zero-mean response assumption were not made, the above relationship would be 
more complicated. Having U now fully defined from equation (20), the modal equivalent of 


equation (17) can be written as 


D D kj\\\,l, ,L (24) 

iDl ,/Dl 

from which can be determined. The matrix on the left hand side of equation (24) involves 4* 

order moments of modal displacements. The right hand side involves 4*, 5* and 6* order modal 
displacement moments since the potential energy is a function of the 2"'^, 3’^'^, and 4^^ order 
displacements. Assuming a Gaussian distributed, zero-mean response, the odd order moments 
are zero and the higher order even moments can be expressed in terms of the 2“'^ order moments. 
For example, 

E § D ^ ^E[q,q, QD E [q.q, \\e ^ .q, § D A [q.q, Qa ^ .q, § . 


Therefore, the matrices of equation (24) can be written solely in terms of the modal response 
covariance. 

4. ITERATIVE SOLUTION 


Because the modal equivalent linear stiffness matrix is a function of the unknown modal 

displacements, the solution takes an iterative form. The time variation of the modal 
displacements and forces may be expressed as: 

q(t) D D qe‘^"‘ f(t) D Q (25) 
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where ('^) indieates the dependency on (o„. Applying (25) to (11) and writing in iterative form 
gives: 

q"" Q (26) 

where m is the iteration number and 

Hmui g QKU[a Kf^ D (27) 

The introduction of the weightings a and P are to aid in the convergence of the solution, with 
the condition that aW pW\. 

For stochastic excitation, (26) is rewritten as: 

5"; D (28) 

and the r,s component of the covariance matrix of modal displacements is 

E[q.q,[S D D S,:„A/7. 

n 

Here, 5"^ is the spectral density matrix of modal displacements and is the spectral density 
matrix of the load in modal coordinates. The diagonal elements, „ , are the variances of the 
modal displacements. is zero for the first iteration, which yields the covariance matrix 
of the linear system. For subsequent iterations {m>\), the nonlinear stiffness depends on 
the minimization approach taken. 


For example, for the force-error minimization approach, for the iteration is determined 
from (14) as 


k: d 



D 


1/7 

r 

I 

I 

1 ^ 


□Z7, [ 

J7 

]D/7i [ 



... 

Ml 


; 








□A 

... E[ 

Ml[ 



(29) 


Substitution of equation (4) into (29) gives 
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( 30 ) 



Since a zero-mean response is assumed, i.e. E[q^ D 0 , equation (30) reduces to: 



(31) 


For the potential energy-error minimization approach, is found by solving the system of 
equations (24) in a similar iterative fashion. 

For both error minimization approaches, the iterations continue until convergence of the modal 
equivalent linear stiffness matrix such that 

The value of s typically used is 0.1%. 

Following convergence, the N x N covariance matrix of the displacements in physical 
coordinates is recovered from 

E . § D D ^ [q^q^ QD ^ (32) 

and root-mean-square (RMS) values are the square roots of the diagonal terms in (32). Further 
post-proeessing to obtain power spectral densities of displacements, stresses, strains, ete., may be 
performed by substituting the converged equivalent stiffness matrix into (11) and solving in the 
usual linear fashion. 

5. IMPLEMENTATION 

The EL proeedures as outlined above were recently implemented within the eontext of 
MSC.NASTRAN [10] version 70.0 (heretofore NASTRAN) using the DMAP programming 
language [15]. Its operation has been verified through NASTRAN version 70.7. Details of these 
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implementations, eolleetively ealled ELSTEP for “Equivalent Einearization using a STiffness 
Evaluation Proeedure,” are doeumented in [16]. The implementations entail first performing a 
normal modes analysis (solution 103) to obtain the modal matrices, from which a subset of L 
modes are chosen. The nonlinear stiffness coefficients are then determined by performing a 
series of linear static (solution 101) and nonlinear static (solution 106) solutions using linear 
combinations of modes as previously described. The iterative solution is performed in a 
standalone routine, which has as its output the RMS displacements in physical coordinates, the 
cross covariance in modal coordinates, and the sum of the linear and equivalent linear modal 
stiffness matrices. The latter can then be substituted for the linear modal stiffness in the modal 
frequency response analysis (solution 1 1 1) for post-processing. 

6. NUMERICAL SIMULATION ANALYSIS 

For validation purposes, a numerical simulation analysis was performed to generate time history 
results from which response statistics could be calculated. The particular method used was finite 
element-based with the integration performed in modal coordinates, as described in [14]. The 
finite element model uses the 4-node, Bogner-Fox-Schmit (BPS) C' conforming rectangular 
element with 24 degrees-of-freedom; 4 bending degrees-of-freedom (w,%,%,°y,,Dy) and 2 
in-plane degrees-of-ffeedom (w,v) at each node. The equations of motion are derived using the 
von Karman large deformation theory. In applying this modal approach, the modal truncation 
should be the same as that used in the EL approach. 

6.1. LOADING TIME HISTORY AND TRANSIENT RESPONSE PROCESSING 

The time history of the load was generated as described in [14]. A time increment (Dt) of 122 
jus was found suitable for the numerical integration and time history records of 2.0 s in duration 

were generated using a record length of 2'“^ samples. A radix-2 number of samples was ehosen 
to faeilitate use of radix-2 EFT algorithms employed for the subsequent analysis. An ensemble 
of time histories was generated by specifying different seeds to the random number generator. 

A typieal time history eorresponding to a pressure load of 8 Pa RMS (1 12 dB) is shown in Figure 
1. The eorresponding probability density function is shown in Figure 2 with the Gaussian 
distribution superimposed upon it. The power spectral density (PSD) for 10 ensemble averages 
gives a speetrum level of 6.25 D lO^^Pa^Hz over a 1024 Hz bandwidth as shown Figure 3. A 
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sharp roll-off of the input spectrum practically eliminates excitation of the structure outside the 
frequency range of interest. 

The structure is assumed to be at rest at the beginning of each loading. An initial transient in the 
structural response is therefore induced before the response becomes fully developed. This 
transient must be eliminated to ensure the proper response statistics are recovered. In the linear 
case, a moving block average of 1.0.S of data was used to ascertain the point in the record at 
which the RMS displacement response became stable. In this manner, it was determined that 
elimination of the first 1.0^ of the 2.0^ record was more than sufficient. While the same criterion 
does not strictly apply to the nonlinear response case because of the dependence of the response 
on the initial conditions, it was employed with satisfactory results. 
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Figure 1; Typical acoustic loading time history. 



Figure 2: Typical probability density of acoustic loading. 
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Figure 3: Typical averaged power spectral density of acoustic loading. 
6.2. RESPONSE STATISTICS 


Response statistics were generated from an ensemble of N=\0 time histories at each load level. 
Estimates of the displacement RMS served as the basis for comparison with the EL method. 
Additionally, confidence intervals for the mean value of the displacement RMS estimate were 
generated to quantify the degree of uncertainty in the estimate [17] using: 


X 


St, 


aDTD 


yfN 


^^n;L/n 

yfN 


ADl 


(33) 


where x and s^ are the sample mean and variance of the RMS estimates from N ensembles, and 
C is the Student t distribution with n degrees of freedom, evaluated at C/2. For the 90% confi- 
dence intervals calculated, C D 0.1 . 


Estimates of the displacement mean, skewness, and kurtosis were also computed to help 
ascertain the degree to which the assumptions made in the development of the EL method were 
followed. Power spectral density and probability density functions (PDF) of the displacement 
were computed for similar purposes. 
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7. RESULTS 


Validation studies were eonducted using a rectangular aluminum plate measuring 0.254m x 
0.3556m of thiekness, h, 0.00102m. All sides were simply supported. The material properties 
used were: 

A D 7.3 D Pa, v D 0.3, p D 2763-^ 
where E is the elastic modulus, v is Poisson’s ratio, and p is the mass density. The plate was 
subjeeted to a spatially-uniform pressure loading over a computational bandwidth of 1024 Hz, as 
shown in Figure 1. Since the loading was uniformly distributed, only symmetric modes were 
included in the analysis. In general, any combination of symmetric and non-symmetric modes 
may be included. 

A NASTRAN model of the full plate was built with 560 CQUAD4 elements measuring 0.0127m 
square for use in the EL analysis. The first two symmetric modes (modes 1 and 4) of this model 
had natural frequencies of 58.38 and 217.27 Hz. For comparison, the natural frequencies given 
by an analytical solution [18] were 58.34 and 216.01 Hz for the first two symmetric modes. 
These two modes were selected as participating modes in the EL and numerical simulation 
analyses. Modal damping was chosen to be sufficiently high (2.0% and 0.54% critical damping) 
so that a good comparison with the numerical simulation results could be made at the peaks of 
the PSD. The finite element model used in the numerical simulation analysis had the same 
element size (0.0127m square) as the NASTRAN model, but a quarter-plate model was used to 
reduce computational time. This model gave natural frequencies of 58.12 and 215.19 Hz for the 
first two symmetric plate modes. Both EL and simulation finite element models were checked 
for eonvergence by running additional analyses with models consisting of 0.00635m elements. 

Analysis was performed at overall sound pressure levels from roughly 106 dB (4 Pa RMS) to 
roughly 160 dB (2048 Pa RMS), in 6 dB increments, giving a dynamic range of 54 dB. Figure 4 
shows the normalized RMS out-of-plane (w) deflection at the plate center as a funetion of 
loading. Both foree-error and potential energy-error minimization implementations are shown. 
The numerieal simulation results are shown with 90% confidence intervals of the RMS estimate. 
At the lowest load level of 106 dB, the response is linear as can be seen by the eomparison with 
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results from a strictly linear analysis (NASTRAN solution 111). Small, but noticeable, 
differences between the linear and nonlinear responses are noted at the 1 1 8 dB load level. The 
degree of nonlinearity increases with load level, as expected. At the highest level of 160 dB, the 
nonlinear response calculations predict RMS center deflections of 2.20 and 2.27 times the 
thickness for the force and potential energy-error minimization approaches, respectively. The 
90% confidence interval on the numerical simulation data by comparison is 2.24 < wrms/^ ^ 2.39 
at 160 dB. Consistent with past observations [11, 12], potential energy-error minimization 
results are closer to known solutions and somewhat higher (by a few percent) than force-error 
minimization results. The linear analysis by comparison predicts center deflections of nearly 13 
times the thickness. 



Figure 4: Normalized RMS center deflection as a function of acoustic load. 

Results from the earlier SEMELRR implementation [9], also shown in Figure 4, were generated 
with the same NASTRAN model. These results are shown to significantly over-predict the 
effect of nonlinearity at loads as low as 124 dB. At 160 dB, the SEMELRR implementation 
predicts an RMS center deflection of 1.38 times the thickness, or approximately 60% less than 
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the mean of the numerieal simulation prediction. The highly non-eonservative nature of the 
SEMELRR analysis makes it unsuitable for nonlinear structural design. 

In order to gain greater insight into the nonlinear dynamics, plots of the time history, PDF, and 
PSD are shown for three load levels, (106, 136 and 160 dB) in the following series of figures. 
Data in the time history and PDF plots correspond solely to numerical simulation results. Data 
in the PSD plots correspond to numerical simulation and EL results, where the EL results were 
generated by running a linear analysis (NASTRAN solution 111) using the modal equivalent 
linear stiffness generated by the EL process previously described. 

Results for the 106 dB excitation level are shown in Figure 5 - Figure 7. This excitation level 
was shown (see Figure 4) to produce a linear response. As expected, the PDF mimics the 
normally distributed PDF of the input shown in Figure 1. The averaged PSD shows excellent 
agreement between the EL, linear, and numerical simulation results. This agreement helps to 
establish the confidence in making comparisons between these two fundamentally different 
analyses. 



Figure 5: Time history of center displacement response at 106 dB. 
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Figure 6; Probability density of center displacement response at 106 dB. 



Figure 7: Power spectral density of center displacement response at 106 dB. 
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Figure 8 - Figure 10 show the nonlinear response associated with the 136 dB excitation level. 
The time history of the center displacement has a visibly higher peak probability and the PDF 
exhibits a minor flattening at the peak. The shapes of the PSDs from the EL analyses are those 
of the equivalent linear systems. Thus, the PSDs from the EL analyses do not show the peak 
broadening effect observed in the numerical simulation. Additionally, harmonics in the 
structural response are present only in the numerical simulation results. The PSDs from both the 
numerical simulation and EL analyses show a positive shift in the frequency of the fundamental 
mode compared with the linear solution. Both EL analyses shift the fundamental frequency by 
nominally the same amount and fall within the broadening fundamental peak of the numerical 
simulation analysis. The frequency of the second mode from the force-error minimization 
analysis also increases and more closely matches the numerical simulation data than does the 
potential energy-error minimization analysis, which shows a negative frequency shift. 



Time (s) 

Figure 8: Time history of center displacement response at 136 dB. 
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Figure 9: Probability density of center displacement response at 136 dB. 



Figure 10: Power spectral density of center displacement response at 136 dB. 
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To investigate this behavior further, the fundamental and second mode frequencies as a function 
of applied load are shown in Figure 11. The shift in the frequency of the fundamental mode is 
consistently positive and close in value between the force and potential energy EL analyses. The 
shift in second mode frequency differs between the two EL analyses, with the force EL results 
showing similar behavior to the fundamental, while the potential energy EL results first 
displaying a drop, then an increase in frequency with increasing load. Since the force-error 
minimization approach fully takes into account the internal force resulting from bending- 
membrane coupling, this approach gives results that are consistent with numerical simulation 
results. The potential energy-error minimization approach places no demand on the internal 
force. Since most of the strain energy in this case is associated with the fundamental mode, that 
behavior is accurately captured. The second mode, having significantly less potential energy, is 
allowed to shift without constraint. 



Figure 1 1 : Shift in fundamental and second mode frequencies as a function of applied load for 

both EL analyses. 


While this behavior appears inconsistent, recall that both EL approaches minimize the error 
between the equivalent linear and nonlinear systems through the modal equivalent stiffness 
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matrix , which is a function of the modal displacement. It is therefore expeeted that only the 

RMS displaeement, or the area under the displacement PSD, should be the same between either 
EL analysis and the numerieal simulation analysis, and not the shape of the PSD itself. This is 
eonsistent with the observations. 


The highest degree of nonlinearity is shown in Figure 12 - Figure 14, corresponding to a 160 dB 
aeoustie load. The time history is further peak oriented and the PDF exhibits substantial 
flattening. The peak broadening in the PSD of the numerical simulation results is severe, and 
nearly flattens the speetrum above 350 Hz. With regard to the EL analyses, positive shifts in the 
frequeneies of the fundamental mode are comparable between the two EL analyses. Shifts in the 
seeond mode frequeneies are more substantial in the force-error minimization approach than in 
the potential energy-error minimization. 
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Figure 12: Time history of center 
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Figure 13: Probability density of center deflection response at 160 dB. 



Frequency (Hz) 

Figure 14: Power spectral density of center deflection response at 160 dB. 
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Moments of the eenter displaeement were calculated from the numerieal simulation results for all 
load levels. They are provided in Table 1 with the RMS center displaeement from both EL 
analyses. The EL and numerical simulation results agree well, thus validating the EL analysis 
over a substantial load range. The validity of assumptions made in the development of the EL 
method is aseertained by observing the mean, skewness and kurtosis. The mean value is 
effeetively zero for all load levels, indicating the assumption of zero mean response has not been 
violated. Although the PDF is more or less skew-symmetric, the shape is flattened at the higher 
load levels as indicated by a decreasing kurtosis from the linear value of 3. The decreasing 
kurtosis values indicate a violation of the Gaussian response assumption. However, even with 
this non-Gaussian response distribution, the EL analyses give good predictions of the RMS 
response. 

Plots of the RMS displacement over the entire plate are shown in Figure 15 - Figure 17 for the 
linear and EL analyses for the 160 dB load case. While the plots appear similar in character, the 
linear displacement contours look rounder than those of the two EL analyses. To better 
understand the nature of these differences, the RMS displacements were normalized with respect 
to their maximum (center) displacement for each of the three cases. The differences between the 
normalized linear and force-error minimization EL displacements, and the normalized linear and 
potential energy-error minimization EL displacements are shown in Figure 18 - Figure 19. 
These are expressed in percent with respect to the maximum normalized displacement (unity). 
In each case, the greatest difference occurs along the horizontal centerline and near the short 
sides of the plate. These plots highlight the need to perform a nonlinear analysis to obtain the 
proper spatial distribution. As previously noted, the difference between the two EL analyses is 
small relative to their displacement, as shown in Figure 20. The difference resembles the seeond 
symmetrie mode. This observation is consistent with the PSD in Figure 14, which indicates the 
maximum magnitude difference between the two EL analyses to be associated with the seeond 
mode. 
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Table 1 ; Moments of the center displacement. 


Load 

(dB) 

Foree Error 
EL RMS (m) 

Potential Energy 
Error EL RMS (m) 

SEMELRR 
RMS (m) 

Num Sim 
Mean (m) 

Numerical Simulation RMS (m) 
(90% Confidence Interval) 

Skewness 

Kurtosis 

106 

2.57xl0'" 

2.57x10'= 

2.56x10'= 

+1.56x10-' 

2.37xlO'=D RMS < 2.76x10'= 

+0.0017 

3.21 

112 

5.12x10'^ 

5.12x10'^ 

5.04x10'^ 

+3 .07x1 O'’ 

d.OSxlO-'^D RMS < 5.47x10'^ 

+0.0013 

3.18 

118 

LOOxlO"* 

1.01x10'“ 

9.54x10'^ 

+1.01x10-'* 

9.34xlO'^D RMS < 1.06x10'“ 

+0.0012 

2.90 

124 

1.88x10'“ 

1.90x10'“ 

1.66x10-“ 

+7.94X10-’ 

1.76xlO'“D RMS < 1.96x10'“ 

+0.0027 

2.60 

130 

3.24x10'“ 

3.30x10'“ 

2.62x10-“ 

+1.87x10-* 

3.42xlO'“D RMS <4.01x10'“ 

-0.0010 

2.81 

136 

5.13x10'“ 

5.27x10'“ 

3.84x10-“ 

-5.23X10-’ 

6.09xl0'“D RMS <6.71x10'“ 

+0.0033 

2.83 

142 

7.68x10'“ 

7.95x10'“ 

5.35x10-“ 

-3.18X10-’ 

8.38xlO'“D RMS <9.30x10'“ 

+0.0032 

2.48 

148 

1.10x10'^ 

1.16x10'^ 

7.27x10-“ 

+2.58x10-* 

1.14x10-^0 RMS < 1.26x10'^ 

-0.0063 

2.38 

154 

1.58x10'^ 

1.63x10'^ 

1.05x10'^ 

+1.04x10-* 

1.61x10-^0 RMS < 1.79x10'^ 

+0.0012 

2.36 

160 

2.25x10'^ 

2.32x10'^ 

1.41x10'^ 

+6.05x10-* 

2.29xlO'^0 RMS <2.44x10'^ 

-0.0061 

2.28 
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1.13E-02 
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5.65E-03 
4.71 E-03 
3.76E-03 
2.82E-03 
1.88E-03 
9.41 E-04 
O.OOE+00 


Figure 15; Contour plot of RMS displacement from linear analysis at 160 dB. 
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Figure 16: Contour plot of RMS displacement from force-error EL analysis at 160 dB. 
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Figure 17: Contour plot of RMS displacement from potential energy-error EL analysis at 160 

dB. 



Percent RMS 
Displacement 
Difference 

5.40 
5.01 
4.62 
4.24 
3.85 
3.46 
3.07 
2.69 
2.30 
1.91 
1.52 
1.13 
0.75 
0.36 
-0.03 



Figure 18: Difference between normalized linear and force-error EL displacements at 160 dB. 
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8. DISCUSSION 


Differences between the two EL analyses and the numerical simulation do exist and warrant 
some discussion. It is seen that the EL approach slightly over-predicts the degree of nonlinearity 
compared to the numerical simulation results. Interestingly, the greatest differences appear in 
the moderately nonlinear regime. The differences do not appear to be due to a violation of the 
assumption of a Gaussian response because the over-prediction does not correlate with 
increasing kurtosis of response. The trends observed here are consistent with comparisons 
between a different numerical simulation analysis and a force-error minimization EL analysis of 
a clamped-clamped beam [13]. Comparisons of the force and potential energy-error 
minimization approaches with an F-P-K solution of a beam, however, indicate that EL solutions 
span the exact solution, with the potential energy-error minimization results being slightly higher 
and the force-error minimization results being slightly lower than the exact solution [12]. It 
would therefore appear that the numerical simulation solutions are less stiff than the EL solutions 
in the moderately nonlinear regime, but this has yet to be fully substantiated. 

Some implications on the use of the EL technique as a basis for fatigue life calculations are 
worth mentioning. First, assuming that stresses or strains from the EL technique will compare 
equally well with those from the numerical simulation analysis, a simple fatigue-life calculation 
based on RMS levels will be much less conservative than calculations based on linear analyses. 
This offers the potential for substantial weight savings for structures designed using nonlinear 
methods. Secondly, it appears that a nonlinear analysis, EL or otherwise, is required to accurately 
ealeulate the RMS deflected shape. Use of a linear RMS deflected shape scaled to the nonlinear 
level would inaccurately reflect the spatial distribution. Simple fatigue-life calculations based on 
the RMS stress or strain could be significantly affected as these quantities depend on the spatial 
distribution of the deformation. Lastly, use of the EL-derived PSD response in a more 
sophistieated fatigue-life calculation requires careful investigation. Reeall that peaks in the 
equivalent linear PSD may oecur at different frequencies than in the PSD from the numerieal 
simulation analysis, as shown in Figure 10 and Figure 14. Methods sueh as speetral fatigue 
analysis [19], whieh take moments of the PSD, may incorrectly account for the eontribution of a 
partieular frequeney eomponent in the cycle counting scheme. It is not known, for example, if 
the narrowly shaped, higher fundamental frequencies of the equivalent linear PSD result in 
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conservative or non-eonservative estimates of fatigue life relative to predietions made using the 
numerieal simulation PSD with more broadly shaped, lower fundamental frequeneies. An 
assessment of this effeet is left as an area for further study. 


9. CONCLUSIONS 

A novel method for determining nonlinear stiffness coefficients of arbitrary struetures has been 
developed. The method can be implemented in any finite element code having a geometrieally 
nonlinear static capability. It has been implemented as a DMAP alter to MSC.NASTRAN to 
demonstrate its effectiveness. 

A potential energy-error minimization EL approach has been extended to handle multiple 
degree-of-freedom systems. The RMS random response predictions from it and the traditional 
force-error minimization approach have been validated through comparison with a numerical 
simulation method over a wide range of load levels. Comparisons with numerical simulation 
results are good, even when the assumption of Gaussian response has been violated. It has been 
shown that a linear analysis grossly over-predicts the RMS displacements (i.e., it is too 
conservative) in comparison with numerical simulation results. The potential energy-error 
minimization approach provided a slightly more conservative estimate of RMS displacement 
than the force-error minimization approach. It was demonstrated that an earlier EL 
implementation (SEMELRR) significantly over-predicts the effect of nonlinearity (i.e., it is non- 
eonservative). 
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